Multiple scattering of light in superdiffusive media 
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Light transport in superdiffusive media of finite size is studied theoretically. The intensity Green's 
function for a slab geometry is found by discretizing the fractional diffusion equation and employing 
the eigenfunction expansion method. Truncated step length distributions and complex boundary 
conditions are considered. The profile of a coherent backscattering cone is calculated in the su- 
perdiffusion approximation. 
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Light transport in disordered media is characterized by 
a multiple scattering process engendered by random fluc- 
tuations of the refractive index in space and described 
perturbatively by expanding the field in powers of the 
scattering potential 1] . Important quantities such as the 
angular distribution of backscattcrcd and transmitted in- 
tensity can be calculated within a diffusion approxima- 
tion and the Green's function of the diffusive equation, 
or diffusion propagator, be used to describe fundamen- 
tal interference effects like speckles and weak localiza- 
tion [H . The diffusive picture has its roots in the concept 
of Brownian motion, where random walkers perform in- 
dependent steps of variable length, each with finite mean 
and variance. By virtue of the central limit theorem, 
the step length distribution after a large number of steps 
approaches a Gaussian distribution irrespectively of the 
microscopic transport mechanism. This rule, however, 
breaks down when the probability to perform arbitrary 
long steps is non- vanishing. In this case, the limiting 
distribution becomes a so-called ct-stable Levy distribu- 
tion 0] and is characteristically heavy-tailed. Such ran- 
dom walks, first studied by Mandelbrot in the framework 
of transport on fractals |3[ , are known as Levy flights if 
steps of arbitrary length can be performed in a unit time 
and Levy walks if performed at a finite velocity As 
transport is dominated by a few huge steps, the resulting 
average distance explored by a walker increases faster 
with time than expected for standard diffusion. This 
type of anomalous transport is called superdiffusion [5( 
and has been found to be ubiquitous in nature [f|. Su- 
perdiffusion of light has recently been observed in het- 
erogeneous dielectric materials [7( and in hot atomic va- 
pors [H . On the theoretical level, superdiffusion has been 
modelled by employing the subordinator method [i| and 
by generalizing the diffusion equation to fractional or- 
der derivatives [l(J. Previous works have evidenced the 



peculiar statistical properties of Levy motions [llj, [12 1 
and shown that several features of real experiments, such 
as properly defined boundary conditions, are nontrivial 
to implement llj], making the description of observable 
quantities nearly impossible. 

In this Letter, we develop a theoretical framework for 
multiple light scattering in superdiffusive media. Our ap- 
proach, which relies on the semi-analytical solution of the 
fractional diffusion equation, allows to study the steady- 



state transport properties of superdiffusive media while 
taking into account the intrinsic finite size of actual ma- 
terials and makes it possible to treat interference effects, 
notably coherent backscattering, in the "superdiffusion 
approximation" . In particular, we calculate the intensity 
Green's function in the superdiffusive regime for vari- 
ous values of the a coefficient and show that arbitrary 
boundary conditions and truncations in the step length 
distribution can be implemented. 

From the microscopic point of view, the disorder aver- 
aged intensity I observed at a point R outside a multiple 
scattering medium can be written as 
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dr i d r 2 dr 3 dr 4 0(r i ) cj)* (r 2 ) ■ 
/(n,r 2 ,r 3 ,r 4 ) G(r 3 , R) G*(R, r 4 ) 
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where cj> is the coherent (i.e. unscattered) propagator 
for the amplitude from outside the sample to the first 
scattering event, G is the averaged propagator for the 
amplitude from the last scattering event to the point of 
observation R (i.e. the solution to the Dyson equation) 
and / (ri, r 2 , r3, r 4 ) is the four-vertex propagator that 
contains all information about transport. When recur- 
rent scattering is neglected, only two two-vertex terms 
contribute to /: the ladder term, which describes the 
incoherent transport, and the most-crossed term, which 
leads to the coherent backscattering cone 13j. More 



complicated combinations of these two terms can also be 
used to describe speckle correlations and intensity fluc- 
tuations 



14(. 



When the step length distribution p{t) decays fast 
enough, the diffusion approximation holds and the prop- 
agator for the incoherent intensity transport satisfies the 
standard diffusion equation [l|. On the other hand, 
for Levy flights, the step length distribution exhibits a 
power-law tail of the formp(^) ~ £-( a+1 ) with < a < 2, 
and the macroscopic motion is described by the diffusion- 
like equation [f|: 



d t C(r,t) = D a V a C(r,t) 



(2) 



where V Q is the symmetric Riesz fractional derivative 
with respect to spatial coordinates and D a is a general- 
ization of the diffusion constant. 

In the normal diffusive case, finite-size effects and in- 
ternal reflection at the boundaries are handled by im- 
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we have: 



FIG. 1. Normalized Green's function / (x,xo) for a slab with 
L = 5, xo = L/3, M = 300, D a = 1 and absorbing bound- 
aries, computed for various values of a. 



posing that the propagator goes to zero at a distance 
from the physical boundaries called the extrapolation 
length 15]. This is possible because the Laplacian opera- 
tor is local in space and thus, the presence of boundaries 
does not change the form of the propagator itself. On 
the other hand, the fractional Laplacian in Eq. [2] is non- 
local and, as such, the superdiffusive propagator is very 
sensitive to the nature of the boundaries. In fact, this 
propagator in non-infinite media is known only in a few 
particular cases 



16|. The description of superdiffusive 



transport in finite-size media therefore requires to deter- 
mine the form of the propagator for arbitrary step length 
distributions and boundary conditions. 

The problem of the non- locality of the fractional Lapla- 
cian can be circumvented by discretizing the fractional 
Laplacian [l7f | , i.e. by replacing the continuous time ran- 
dom walk by discrete hops on a lattice and V" by an 
M x M matrix, which, when applied to the vector repre- 
senting our function C, converges to the continuum oper- 
ator when M goes to infinity. Let us consider a ID system 
where space and time are discretized (the generalization 
to the 3D case will be presented later in this Letter) and 
define as the probability to perform a jump from 

site i to site j. The macroscopic transport is expected to 
be described by Eq.[2]when the time and space discretiza- 
tions are fine enough. Concurrently, the microscopic re- 
distribution process that occurs at each time interval can 
be written as Cj(t n+ i) = X)i ^Ij-ilQ^n) such that [l^ |: 



Cj(t n+ \) - Cj(t n ) 



1 M 

~5> 



5 t , 3 )C,j{t n ), (3) 



where Si^ is the Kronecker delta. 

In the limit r — > 0, the left- and right-hand sides of 
this equation represent dtCj(t) and D a V a Cj(t), respec- 
tively. Any microscopic redistribution property, provided 
that it leads to superdiffusion, can therefore be used to 
write a discretized version of the fractional Laplacian. In 
particular, if il is the matrix of transition probabilities, 
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(4) 



The convergence properties of this limit depend on the 
particular choice of the transition probabilities. The 
most natural choice, Ct?u_ji oc \i — j\ ' Q+1 ^ is known to 
suffer from a poor convergence, especially when a — » 
2 A much faster convergence has been demon- 

strated [13] using a direct discretization of the fractional 
Laplacian, leading to 17]: 
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where h is the distance between two consecutive nodes 
on the lattice and F is the Euler gamma function. In 
this framework, setting = when j is outside a 

given interval [0, L] corresponds to the situation in which 
walkers stepping out of the interval cannot ever re-enter 
it. Thus, reducing the infinite size matrix fi to a M x M 
matrix comes to imposing a finite size with absorbing 
boundary conditions to the system (20j . 

A physical model for superdiffusion in real systems 
should rely on Levy walks rather than on Levy flights 
since all jumps are bound to have a finite velocity. The 
resulting spatiotemporal coupling 4] is known to make 
the description of Levy walks difficult to handle analyti- 
cally [2l| as opposed to Levy flights, essentially described 
by Eq. [2] This coupling, however, becomes irrelevant in 
the steady-state regime since the amount of time required 
to perform a jump plays no role. The intensity Green's 
function in a ID superdiffusive medium for a continuous 
point source at xq is then given by the following time- 
independent fractional differential equation: 



D a y a f (x, x Q ) = -5{x - .t ). 



(6) 



Note that f (x, Xq) is the two- vertex propagator appear- 
ing in Eq. [TJ in the superdiffusion approximation. 

A complete description of the operator V" is provided 
by the eigenfunctions ipi and eigenvalues of the matrix 
fl — 1, with h = L/AI. Then, / can be expanded in 
terms of the eigenfunctions ipi as a linear combination 
f ( x ) — ^2i a i' l Pi) where ais are the coefficients to be 
determined. By substituting the above expressions in 
Eq. [S] and considering that ipj(x)'ijji(x) = Si.j, we obtain 
ai = —ipi(xo) (D a Xi) , yielding the Green's function of 
a ID superdiffusive motion of exponent a: 



M 

E 



tjji (x) Ipj (Xq) 

D a Xi 



(7) 



Figure [T] shows the normalized Green's function / for 
different values of a. For a = 2 we recover the familiar 
triangular shape typical of the diffusive regime while the 
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FIG. 2. Normalized Green's function / (x,xo) for a slab with 
a = 1.5, L = 5, xo — L/3, M = 300 and D a — 1, computed 
for various values of (a) i max and (b) the reflectance R on the 
left boundary. 



Green's function becomes more and more cusped when 
a decreases. 

The microscopic transport mechanism of light in real 
systems is subject to additional and more complex fea- 
tures, including a truncation in the step length distribu- 
tion and partially or totally reflecting boundaries. Those 
can be taken into account through Eq. [4] by modifying 
the transition probabilities in the medium. Trun- 

cated Levy step distributions [22| , which are unavoidable 
in finite-size systems jjj], can be implemented by setting 
to zero all transition matrix elements where \i — j\ > l max 
and renormalizing fi so that = IVi before in- 

troducing the boundaries. Figure [2^, shows how the nor- 
malized Green's function at constant a changes with the 
truncation length. While when / max ~ L there is only 
a minor correction to the shape of the Green's function, 
f(x,Xo) becomes very similar to the diffusive one when 
L 



< L. 



Partially reflecting boundaries are also expected to 
change the shape of the propagator. They can be imple- 
mented by considering that walkers reaching a boundary 
of the superdiffusive media have a probability R to be 
reflected. In the case where the left boundary has re- 
flectance R while the right one is totally absorbing, par- 
tial reflection is enforced by mapping all matrix elements 
corresponding to j < onto their mirror images, yielding 



n 



1,3 



Rl 



The case of both partially 



reflecting boundaries are conceptually analogous albeit a 
bit more involved due to the fact that one has to con- 
sider the possibility of performing very long steps that 



might bounce forth and back a large number of times. 
Figure (2b shows the normalized Green's function for dif- 
ferent values of the reflectance R of the left boundary. 
Note that when R = 1 the gradient of f(x,Xo) on the 
left boundary goes to zero as expected [16| . 

Up to now, we have considered the ID case of superdif- 
fusive transport in finite media. Our approach can easily 
be extended to higher dimensions in the case of a slab 
geometry. Orienting the slab such that its interface is 
normal to the re-axis, the system becomes translationally 
invariant in both y- and z-directions. The 3D counter- 
part of Eq. [6] can then be written in terms of the Fourier 
transform of /(r,ro) in the yz-plane as: 

D a (VS - kl) f (x, x , k x ) = -5(x - x ). (8) 

Applying to Eq. [8] the approach used to find Eq. [71 we 
find: 



f(x,x ,kx) = -^2 



ipi(x )ipi{x) 
i D a (Xi — kjj 



(9) 



where k± = |kj_|. The 3D Green's function of the su- 
perdiffusive medium can be obtained at this point by 
performing an inverse Fourier transform of Eq. [9] The 
full intensity distribution in the system (including the 
transmission profile .7]) can be obtained upon integration 
over a suitable source. As shown below, Eq.|H]can also be 
used to compute the shape of the coherent backscattering 
cone in the superdiffusion approximation. 

In the multiple scattering regime, interferences can 
play a major role in determining the transport properties 
of the medium. In particular, in the exact backscattered 
direction, the interference coming from counterpropagat- 
ing paths is always constructive as long as the system is 
reciprocal [lj. This leads to a narrow cone of enhanced 
albedo in reflection known as the coherent backscatter- 
ing cone. If reciprocity is not broken, the peak of the 
enhanced albedo is exactly twice the common incoher- 
ent reflection and presents a triangular cusp on the top 
(23| while its exact shape depends on transport in the 
medium. As a rule of thumb, long paths contribute to 
the formation of the cusp and short ones to the tails of the 
cone. Since in a superdiffusive regime there is no a priori 
reason for reciprocity to break down, we expect longer 
paths to contribute more than in standard diffusion and 
thus, expect a sharper peak. 

The coherent component of the albedo A can be calcu- 
lated starting from Eq. Q] by setting ri = r4 and r-2 = r$ 
[lj. Considering a planewave at normal incidence on the 
slab interface, using the Fraunhofcr approximation for 
the Green's functions and assuming that the step distri- 
bution follows a power-law, we can write: 



( ri ) = ^ (a+1)/2 ^ 
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FIG. 3. Normalized coherent albedo for a slab with L = 5, 
M — 300 and D a — 1, computed for various values of a. 
When a is decreased the contribution of long paths to the 
backscattering cone is larger leading to a sharper cusp and 
to tails that fall down more slowly (but still go to zero for 
8 — ±7r/4). A zoomed view of the top of the curves is shown 
in the inset. 



and similarly for </>*(r2) and G*(R, ri), where 8 is the 
angle with respect to the normal of the slab, kj and k c 
the wavevectors of the incident and emergent planewaves 
respectively, and kj^ = (k; + k e )j_. After substitution in 
Eq. [I] and Fourier transform, we obtain: 



A oc JJ dxidx 2 (~|) f(xi,X2,kx)- 



(10) 



Using the discretized approximation in Eq.[5]for the prop- 
agator f(xi,X2, kj_), we find the following expression for 
the coherent albedo from a superdiffusive slab: 



A oc 
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M 



^D a (\ t -\ksm8\ a )- 



(11) 

Figure [3] shows the normalized profile of the coher- 
ent backscattering cone as a function of a. When a 
is decreased the amount of light transmitted through 
the sample is increased (in the superdiffusive regime 
T oc L~ a l' 1 [l9j]) yet, at the same time, long steps be- 
come increasingly important, thereby making the profile 
more cusped. We note that a similar effect has been pre- 
dicted for enhanced backscattering in fractal media 24 1 . 
Note also that since all calculations are done considering 
a finite thickness, part of the light is lost by transmission 
through the system, and thus, the top of the cone for 
a = 2 appears rounded 

In conclusion, we obtained a semi-analytical formula- 
tion for the intensity Green's funtion of multiple scattered 
light in the superdiffusive regime applying the eigenfunc- 
tion expansion method to the discretized version of the 
steady-state fractional diffusion equation. This approach 
makes it possible to describe the behavior of many ob- 
servable properties of superdiffusive media of finite size 



with complex boundary conditions (absorbing, partially 
reflecting, reflecting) as well as truncated step distribu- 
tions. It also allows for the calculation of fundamental 
interference effects, such as the coherent backscattering 
cone, in the superdiffusion approximation. 
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